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Abstract. It is customary to perform analy- 
sis of the Earth's rotation in two steps: first, 
to present results of estimation of the Earth ori- 
entation parameters in the form of time series 
based on a simplified model of variations of the 
Earth's rotation for a short period of time, and 
then to process this time series of adjustments 
by applying smoothing, re-sampling and other 
numerical algorithms. Although this approach 
saves computational time, it suffers from self- 
inconsistency: total Earth orientation parame- 
ters depend on a subjective choice of the apriori 
Earth orientation model, cross-correlations be- 
tween points of time series are lost, and results of 
an operational analysis per se have a limited use 
for end users. An alternative approach of direct 
estimation of the coefficients of expansion of Eu- 
ler angle perturbations into basis functions is de- 
veloped. These coefficients describe the Earth's 
rotation over entire period of observations and 
are evaluated simultaneously with station posi- 
tions, source coordinates and other parameters 
in a single LSQ solution. In the framework of 
this approach considerably larger errors in apri- 
ori EOP model are tolerated. This approach 
gives a significant conceptual simplification of 
representation of the Earth's rotation. 
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for the Earth's rotation is not only a subject 
of academic interest, but is the area for vari- 
ous important applications. The complexity of 
the formalism increases the probability of misun- 
derstanding. Misunderstanding significantly in- 
creases the risk of an error in a software which 
implements reductions. In the worst case sce- 
nario, for instance, a glitch in a navigation equip- 
ment of a passenger aircraft controlled by such 
a program may result in a wreck and loss of life. 
That is why it is vitally important to simplify 
description of the Earth's rotation for end users. 

An alternative approach for describ i ng th e 
Earth rotation was proposed bv IPetrovl ( 20061 ). 
It was demonstrated that instead of using a very 
complex mathematical model for the a priori ro- 
tation matrix and time series of the Earth ori- 
entation parameters which correct that a priori 
rotation matrix, it is feasible to represent the 
Earth's rotation in a form of a sum of a very 
simple, coarse a priori mathematical model and 
a set of coefficients of expansion of the pertur- 
bational rotation into the Fourier and B-splinc 
bases which are estimated from observations. 

In this paper the implications of the alterna- 
tive approach with respect to the traditional ap- 
proach, its advantages and limitations, are dis- 
cussed. 

2 Definition of the Earth's rotation 



1 Introduction 

The complexity of the traditional formalism 
for describing the Earth's rotation is frustrat- 
ing. Even in the community of professional as- 
tronomers a complete understanding of the pro- 
cedure of applying reduction for the Earth's ro- 
tation in all details is not common, especially of 
the last IAU 2000 IAU 2006 recommendations 
(jMcCarthv &; Petit] . |2004) . However, reduction 



Space geodesy techniques, such as GPS, SLR, 
DORIS and VLBI allow to measure the time 
intervals or the differences of time intervals of 
electromagnetic wave propagation from observed 
bodies to observing stations. Solutions of dif- 
ferential equations of wave propagation ties po- 
sition vectors of observed bodies and their time 
derivatives with position vector of observing bod- 
ies and their time derivatives. Therefore, a rela- 
tive motion of observing stations with respect to 
observed bodies can be determined from analy- 



sis of observations. Since the observing stations 
are fixed with respect to the Earth's crust, it 
is convenient to express their positions in a ter- 
restrial coordinate system in which a motion of 
stations would be small, an order of magnitude 
of Earth's crust deformation, i.e. 1CP 8 . A time- 
independent reference position vector in such a 
coordinate system, r T is determined via the fol- 
lowing linear transformation: 

r(t) c = M a (t) r T + q e (t) X r. T + d a (t) T + d e (t) T (1) 

here f(t) c is the reference position vector in ce- 
lestial coordinate system, d a {i)T and d e (t)T are 
the time-dependent a priori and empirical station 
displacement vectors, M a (t) is the a priori rota- 
tion matrix, q e (t) is the empirical small vector 
of perturbational rotation. It should be noted 
that such a transformation does not determine 
coordinates of position vectors in the terrestrial 
coordinate system uniquely. Equation [1] should 
be accompanied with the vector equation of ad- 
ditional boundary conditions: 
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q e (t) x d e (t)i — const 



(2) 



where summing is performed over all stations of 
the network. The rotation matrix A4 a (t) and its 
perturbation q e (t) can be considered as the ma- 
trix of Earth's rotation if a) the relative motion 
of stations with respect to the local areas of the 
Earth's crust is negligible; b) the reference posi- 
tions of observed bodies are at reset in an inertial 
coordinate system. Thus, the Earth's rotation is 
defined here as a net rotation of a polyhedron 
of reference positions of observing stations with 
respect to a polyhedron of reference positions of 
observed bodies. 

This definition of the Earth's rotation follows 
the materialistic approach: "the Earth's rota- 
tion is that what is observed". The currently 
prevailing paradigm tends to consider this phe- 
nomena from the idealistic point of view, which 
involves considering a relative motion of celes- 
tial spheres, one of which represents the ideal 
Earth and an other represents the sky (refer to 
(lMoritj . [l987h more details) That approach op- 
erates such notions as big circles, poles, axes, 
ecliptic, true equator of date etc., and follows 
the tradition which can b e tra ced back to an- 
cient Greece (|PannekoekL Il989h . This concept 
has a mechanical inter pretation of a rolling and 
sliding Poinsot's cones ( PoinsotJ . 1834). However, 



this concept is not adequate when one needs to 
describe the Earth's rotation with an accuracy 
comparable to the precision of modern obser- 
vations, i.e. 3-5 • 10~ 10 rad. First, the Earth 
can not be considered as a rigid body, and a 
notion of the rotation axis, i.e. a set of points 
to which a distance from any point that belongs 
to the body does not change with time, is not 
applicable any more. Second, the presence of 
high frequency variations in the Earth's rotation, 
others than retrograde diurnal terms, makes the 
Poinsot's interpretation inadequate. Thirdly, the 
complexity of the mathem atical model — t he so- 
call MHB2000 expansion (|Mathews et all l2002h 
consists of more than one thousand separate mo- 
tions — makes a meaningful interpretation prob- 
lematic, similar to the situation of description 
of the planetary motion in terms of Ptolemy's 
epicycles. 

3 Parameterization of the Earth's 
rotation 

Since both M (t) and vector d(t) are functions of 
time, i.e. infinite sets of points, their evaluation 
from a finite set of observations can be performed 
only in the form of representing them via some 
functions. The choice of these functions we call 
a mathematical model. The model depends on a 
finite set of unknown parameters which are de- 
termined from observations. 

The fundamental problem is that no observa- 
tion technique, except the laser gyroscope, is sen- 
sitive to the instantaneous Earth's rotation vec- 
tor or its time derivatives directly. The rotation 
angles can be derived using the least square es- 
timation procedure together with evaluation of 
other parameters, and it requires the accumula- 
tion of the sufficient amount of data in order to 
separate variables. The estimates of the Earth's 
rotation angles cannot be sampled too fast. A 
typical sampling rate of estimates is one day, 
since this allows compensation of a certain type 
of systematic errors. In some cases the sampling 
rate can be reduced to several hours. 

3.1 Traditional approach 

A careful examination of th e traditional 
appro ach first proposed by Irlerring et al. 

1986), reveals th ree mathematical models 

McCarthy fc Fetid . |2004 ): 1) the a priori 



model, 2) the estimation model, and 3) the post- 
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processing model. Corrections to rotation angles 
around coordinate axes Ai(t) are parameterized 
in the form equivalent to 

Ai(t) = c(t) cos— Cl n t + s(t) sin— Q, n t + 

bx{t) + b^^it-to) 

A 2 (t) = c(t) s\n-Q, n t - s(t) cos -VL n t + (3) 

b 2 (t) + b 2 (t)*(t-t ) 

A 3 (t) = bait) + b 3 (t)*(t-t ) 

where Q n is the Earth's angular velocity. 

This estimation model is not adequate for 
a long period of time. Usually parameters 
s,c,b,b are determined for a given 24 hour pe- 
riod. Then the time series of Sj, q, bi are filtered 
and smoothed. The coefficients of the interpolat- 
ing spline are computed using the smoothed time 
series of these parameters. This interpolating 
spline for continuous functions s(t),c(t),b(t) de- 
fined at the entire interval of observations forms 
the post-processing model. The a priori rota- 
tion matrix -M a (t) is represented as a product 
of 12 elementary rotations around coordinate 
axes. Some angles of these rotations depend on 
functions s(t), c(t),b(t) from the post-processing 
model. 

The traditional approach has certain disad- 
vantages. The three mathematical models in- 
volved contradict each other. The parameter es- 
timation is not optimal: the raw time series min- 
imizes residuals in the least square sense, but the 
filtered and smoothed series does not; if it did, no 
smoothing would have been needed. The estima- 
tion modcl[3]has a very small range of validity. In 
order to compensate its coarseness, a very refined 
and sophisticated a priori model is required. It 
should be accurate at a very high level. Small 
changes in the a priori model results in changes 
not only adjustments, but in total Earth orienta- 
tion parameters due to inadequacy of the estima- 
tion modelEl To demonstrate it, a series of VLBI 
solutions was made. The first reference solution, 
gsf2006d3, used as a priori the USNO Final time 
series of pole coordinates and UT1 with the time 
span of 1 dajH. Eight parameters, s, c, b, b, were 
estimated for each 24 hour observing session in- 
dependently. In the second solution the Gaus- 
sian noise with the standard deviation 1 nrad was 
added to the USNO Finals EOP. In other trial 
solutions the Gaussian noise was passed through 
the rectangular low-pass digital filter with the 

J http : //vlbi .gsfc .nasa.gov/soltions/2006d 
2 f tp: //maia.usno .navy .mil/ser7/f inals . all 



Table 1. The rms of the differences in total EOP of 
solutions with a priori time series with added Gaus- 
sian noise with respect to the reference solution in 
nrad. The Gaussian noise was passed though the 
rectangular filter which cut the frequencies that cor- 
responding periods shorter than a threshold. The 
rms of the noise was 1 nrad in all cases. 



Freq. cutoff 


Pole coord. 


UT1 


Nutation Ae 


Periods < l d 


0.147 


0.138 


0.148 


Periods < 3 d 


0.082 


0.074 


0.084 


Periods < 7 d 


0.048 


0.021 


0.016 


Periods < 10 d 


0.041 


0.021 


0.008 


Periods < 15 d 


0.038 


0.012 


0.004 



frequency cutoffs which correspond to periods of 
3, 7, 10 and 15 days. The noise was re-scaled in 
order to have the standard deviation of 1 nrad in 
all cases. The rms of the differences in totals, the 
sum of apriori and adjustments, with respect to 
the reference solution are presented in table [1] 

This factor is one of the reasons of the so- 
called "analysis noise" — the differences in re- 
sults of analysis of the same data processed by 
different analysis centers. In order to reduce 
the influence of errors in a priori to a negligible 
level, say 1/2 of the formal uncertainty, which is 
nowadays is at a level of 0.3 nrad, the errors in 
the a priori EOPs must be very small, less than 
1 nrad. 

Another disadvantage of the traditional ap- 
proach is that results of parameter estimation are 
not usable directly by an end-user: they should 
be post-processed. This makes it difficult to as- 
sess the errors of the interpolated time series of 
c(t), s(t), b(t). Correlations are lost, contribution 
of errors of the a priori EOP model is not taken 
into account. This makes evaluation of statisti- 
cal hypothesis based on estimates of the Earth 
orientation parameters problematic. 

3.2 Alternative approach: expansion of the 
perturbational rotation into basis 
functions 

We can overcome the setbacks of the traditional 
approach if we refine the estimation model and 
make it valid not only over a 24 hour period, but 
over the entire time range of observations, i.e. 
several decades. This can be done in a form of ex- 
pansion of the vector of perturbational rotation 
q e (t) into series over basis functions. The choice 
of basis functions is determined by the nature of 
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the process under investigation. It is proposed 
to expand q e (t) into Fourier and B-spline bases 
in this paper: 
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where B™ (t) is the B-spline function 

( Niir nb er geit 1 1 9 8 9h of degree m determined at a 
sequence of knots k-m, • • • , h, h, ... t k ; 

ujj are the frequencies of external forces; 
the coefficients f lkl P^, P^ S c , S s , E^ E^ are 
the parameters of the expansion; Cl n is the 
nominal frequency of the Earth's rotation 
7.292 115 146 707 • 1CT 5 rad s _1 . 

The B-spline basis is adequate for modeling 
a smooth, slow component in q e (t), the Fourier 
basis is adequate for modeling harmonic varia- 
tions, the coefficients of the cross-term, S C ,S S 
take into account corrections to the precession 
rate and the ecliptic obliquity rate. 

Although the B-spline and Fourier bases alone 
are orthogonal, the sum of two bases is in general 
not orthogonal. The estimation model should be 
complemented by the orthogonality condition: 
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The choice of the degree of B-spline and the 
time span between knots, the number of con- 
stituents and the frequencies of the Fourier basis 



elements are determined by the targeted accu- 
racy of the estima t ion m odel. It was demon- 
strated by iPetrovi (l2006h . that the time span 
3 days for qi (t) , q% (t) , 1 day for the qs(t) com- 
ponent, and 740 harmonic constituents is suffi- 
cient for keeping the accuracy of the estimation 
model below a lO^ 11 rad level. The frequencies 
of the Fourier constituents are selected in such a 
way that they correspond to excitation caused by 
a) tide-generating potential; b) the near diurnal 
free wobble; c) the atmospheric nutations. 

There are several approaches for selecting the 
frequencies of relevant Fourier constituents. In 
the first approach a theory of the Earth's nuta- 
tions is used. Since the spectrum of the tide gen- 
crating potential consists of the series of sharp 
peaks, including constituents with the ampli- 
tudes of the potential and forced nutations higher 
than a certain limit is sufficient. The near di- 
urnal free wobble and the atmospheric nutation 
are band-limited processes, so selecting all the 
frequencies within those bands with a step recip- 
rocal to the length of observations is sufficient 
for modeling these components. 

Some components of the theoretical spectrum 
have a frequency separation less than 1/T, where 
T is the time interval of observations. In this 
case strong constraints are imposed to force the 
ratio of complex amplitudes to be the same as 
theoretical. 

Though exploiting theoretical knowledge 
about the Earth rotation allows to reduce consid- 
erably the number of constituents of the Fourier 
basis, this raises a certain concern. If the estima- 
tion model implicitly incorporates theoretical as- 
sumptions, strictly speaking these estimates can- 
not be used for validation of the theory. Al- 
though the ratios of amplitudes between close 
constituents of EOP variations caused by the 
tidal potential exerted by external bodies have a 
sound theoretical basis, we should bear in mind 
that the ultimate goal of comparison of theoret- 
ical predictions with observations is to check va- 
lidity of assumptions put in the foundation of 
the theory and to make a judgment whether the 
model is complete or not. If there are unac- 
counted additive constituents at these frequen- 
cies, for example, caused by the free motion, by 
the atmospheric or oceanic excitation, the theo- 
retical ratios may not be valid. 

Anther way to find the sequence of frequencies 
where the signal is present is to make a set of trial 
solutions and to estimate amplitudes of harmonic 
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variations of q e (t) at a set of frequencies sam- 
pled with a step reciprocal to the length of ob- 
servations within the diurnal, semi-diurnal and 
ter-diurnal bands and discard the constituents 
where no statistically significant signal was de- 
tected. This approach is free from a potential 
bias of adopted theories. 

4 Estimation of the vector of the 
perturbational rotation from VLBI 
data 

The VLBI dataset from January 1984 through 
September 2006 was used for validation of the 
proposed approach. On average, 150 twenty 
four experiments per year have been observed. 
The number of participating stations in each in- 
dividual session varies from 2 to 20, although 4-7 
is a typical number. No station participated in 
all sessions, but every station participated in ses- 
sions with many different networks. All networks 
have common nodes and, therefore, are tied to- 
gether. 

The requirement to the accuracy of the a pri- 
ori Earth's rotation model is determined by a 
condition that ||<f e || 2 be less than errors of de- 
termination of the perturbational vector q e , i.e. 
10~ n rad. This gives the requirement of the ac- 
curacy of the a priori matrix: 3 • 10 rad. This 
is three orders of magnitude lower than the re- 
quirement to accuracy of the a priori rotation 
matrix in accordance with the traditional ap- 
proach. We can exploit this advantage of the al- 
ternative approach and use the simplest possible 
model. The following expression for the a priori 
matrix of the Earth's rotation A4 a (t) according 
to the Newcomb-Andoyer formalism was used: 



M a (t) =ft 3 (Co) • K 2 (9o) ■ U 3 (z) ■ K^-eo) 
TZ 3 (A^) ■ n x {e a + Ae) • U 3 (-S) 



(6) 



where TZ i is a rotation matrix around the axis 
i. For the variables Co, 9q, z , e o, Aip, Aeo, 5*, the 
following simplified expression were used: 



Co = Coo + Coi t + Co2 t 
#0 = #00 + #01 t + #02 t 2 

z = z a + z\ t + z 2 t 2 

e = e 00 + C01 t + £02 t 2 
2 

Aip = pj sin (a.j + (3j t) 



Ae = H e J COS (°3 + #/ l ) 

3 

S = So + E Q + Aip cos e 

+ +Coi +2i +E{)t 

+ (C02 +Z 2 + E 2 ) t 2 
2 

+ Y ( E * cos ^ f + E t sin 7* t) 



Here t is TAI time elapsed from 2000, January 
01, 12 hours. Some of these parameters were 
taken from theory, some of them were found with 
LSQ fit of time series of adjustments of pole co- 
ordinates and UT1. Errors of this apriori Earth's 
rotation model are less than 2 • 10 -6 rad over the 
period 1984-2006. 

Unlike to processing GPS observations, anal- 
ysis of VLBI observations is done in a so-called 
global mode: a set of 2,000-20,000 global pa- 
rameters which are considered common over the 
entire set of observations, i.e. 22 years; 1,000- 
50,000 local parameters specific for a given 24 
hour session; and over a million of segmented 
nuisance parameters specific for a 20-60 minute 
interval are solved in a single least square solu- 
tion directly. 

In the present solution global parameters were 
stations positions, station velocities, amplitudes 
of harmonic variation in site positions, coeffi- 
cients of B-spline modeling non-linear site po- 
sition variations, source coordinates and proper 
motions, harmonic variations at 740 frequencies 
and coefficients of B-spline of the perturbational 
vector of the Earth's rotation q e . 

4.1 Summary of the VLBI results 

Results of analysis are available on the Web 
at http://vlbi.gsfc.nasa.gov/erm The 
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weighted root mean square of postfit residuals 
was the same as in the solution that followed the 
traditional approach, 21.9 ps. 

For comparison with the USNO time series of 
pole coordinates and UT1-TAI, these time se- 
ries were transformed to the vectors of perturba- 
tional rotation, and the coefficients of the cu- 
bic interpolating spline was computed. These 
coefficients represents the USNO mathematical 
model of the Earth rotation. The accuracy of 
the empirical Earth rotation model is higher at 
the instants of time of middle of observing ses- 
sions. The values of vectors of perturbational ro- 
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Table 2. The first row shows the rms of the dif- 
ferences over the period 1996.0-2006.0 between two 
models of the Earth's rotation: the empirical model 
and the USNO model. The second row shows the rms 
of the differences between solution gsf2006d that fol- 
lows the traditional approach and the USNO model. 



Solution 


Pole 


UT1 


Pole rate 


UT1 rate 




ht 9 


rad 


10 -14 


rad s _1 


ERM 


0.64 


0.52 


0.99 


0.81 


gsf2006d 


0.56 


0.42 


1.95 


1.52 



tation and their time derivatives were computed 
at these moments of time for both the empiri- 
cal model from the VLBI solution and from the 
USNO mathematical model. The rms of these 
differences are presented in table [2] It should be 
noted that the GPS estimates of pole coordinate 
almost entirely dominate the USNO time series. 
Therefore, the differences between pole coordi- 
nates and their rates give us a measure of an 
agreement between the empirical Earth rotation 
model from VLBI and the independent estimates 
from GPS. 

In order to validate the estimates of the har- 
monic variations of the perturbational rotation 
vector, a trial solution following the traditional 
Earth rotation parameterization and the a pri- 
ori empirical model of the Earth rotation from 
the previous solution was run. The rms of ad- 
justments of nutation angles over 1996.0-2006.0 
with respect to the aprio ri MHB2000 expan- 
sion ([Mathews et all I2002T ) was 0.98 nrad and 
0.39 nrad with respect to the empirical Earth 
rotation model. 



5 Conclusions 

It was demonstrated that the empirical Earth 
rotation model can be determined directly from 
observations over a period of 22 years using the 
least square estimation technique. The advan- 
tage of the proposed approach is that a con- 
tinuous function describing the Earth's orien- 
tation is determined in one step without pro- 
ducing intermediate time series. Another ad- 
vantage of the proposed approach is that a 
simplified a priori model with only 31 numeri- 
cal parameters is sufficient, while according to 
the traditional approach a complicated a pri- 
ori model of precession, nutation, high frequency 
harmonic variations of the Earth's rotation and 
a filtered and smoothed time series of the Earth 



orientation parameters produced in the previous 
analysis, in total 46 000 numerical parameters 
( McCarthy fc Petitl . Eiol . is needed. 

The proposed approach is conceptually much 
more simple than the traditional approach, since 
it does not operate with idealistic notions that 
are not observable, such as the non-rotating 
origin, the equinox, various intermediate poles, 
axes, etc. 

The EOP, the station position and velocities, 
the source coordinates are produced in a sin- 
gle LSQ adjustment, and therefore, are mutually 
consistent. 

It was demonstrated that the empirical Earth 
rotation model derived from analysis of VLBI ob- 
servations gives the differences with respect to 
the EOP derived from analysis of independent 
GPS observations at moments of observation at 
the same level, within 20%, as the differences of 
the VLBI EOP series produced with the tradi- 
tional approach. The advantage of the proposed 
approach is that the estimates of the EOP rates 
are by a factor of 1.5-2.0 closer to the GPS time 
series than the VLBI EOP rates estimated fol- 
lowing the traditional approach. The rms of esti- 
mates of residual nutation angles with respect to 
the empirical Earth rotation model is 2.6 times 
smaller than the residuals with respect to the 
MHB2000 expansion. 
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